Quantum effects in the dynamical localization of Bose-Einstein condensates 

in optical lattices 
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We study quantum effects in the dynamics of a Bose-Einstein condensate loaded onto the edge of 
a Brillouin zone of a one-dimensional periodic potential created by an optical lattice. We show that 
quantum fluctuations trigger the dynamical instability of the Bloch states of the condensate and 
can lead to the generation of arrays of matter-wave gap solitons. Our approach also allows us to 
study the instability-induced anomalous heating of the condensate at the edge of the Brillouin zone 
and growth of the uncondensed atomic fraction. We demonstrate that there are regimes in which 
the heating effects do not suppress the formation of the localised states. We show that a phase 
imprinting technique can ensure the formation of gap soliton trains after short evolution times and 
at fixed positions. 

PACS numbers: 03.75.Lm 



I. INTRODUCTION 

Bose-Einstein condensates (BECs) loaded into optical 
lattices provide a unique opportunity for testing many of 
the fundamental concepts of solid state physics in a flex- 
ible and defect-free model system. However, one feature 
sets this system apart from any condensed matter system 

the intrinsic nonlincarity of the BEC that is fundamen- 
tally due to elastic atomic scattering. The physics of the 
intricate interplay between the periodicity of the lattice 
and nonlinearity of coherent matter wave has recently 
been subject to numerous theoretical and experimental 
studies, see Refs. 0, Q for an overview. 

One of the most striking manifestations of this inter- 
play observed so far is formation of bright atomic solitons 
in a condensate with repulsive nonlinear interactions 0. 
This unique form of nonlinear localisation occurs within 
the gaps of the periodicity-induced band-gap spectrum 
of the matter Bloch waves, and is possible only due to 
anomalous diffraction properties that the matter waves 
acquire at the edges of the Brillouin zone due to the Bragg 
scattering on a periodic potential. The experimental ob- 
servation of a localized excitation (gap soliton) required 
preparation of the condensate with a small number of 
atoms in the ground state of the one-dimensional lattice 
in the middle of the Brillouin zone. The BEC was then 
adiabatically driven to the band-edge in an accelerating 
lattice, followed by evolution at the band edge in the 
lattice moving with a constant velocity. 

Employing a similar procedure, but with condensates 
containing a large number of atoms, Fallani et al. Q ob- 
served another basic nonlinear effect — the dynamical 
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(modulational) instability of a Bloch wave evolving at 
the edge of a Brillouin zone in a moving lattice. The ob- 
served signatures of this instability were significant loss of 
atoms from the condensate and its spatial fragmentation. 
The loss of atoms is manifested in the "anomalous heat- 
ing" , i.e. enhanced growth of the thermal or uncondensed 
fraction of atoms. On the other hand, theoretical studies 
of condensates with a large number of atoms at the edge 
of the Brillouin zone beyond the onset of the dynami- 
cal instabilities have shown that the instability-induced 
dynamics can also result in condensate localization and 
formation of a train of the localized gap soliton-like struc- 
tures p| . These studies were performed in the framework 
of the mean-field (Gross-Pitaevskii) model, and therefore 
the loss of atoms from the condensate and formation of 
a thermal cloud could not be captured by this analysis. 



The purpose of this paper is to answer the ques- 
tion: Can periodic condensate localization at the band 
edge still occur in the presence of anomalous heating? 
To this end, we perform an analysis of the condensate 
non-adiabatically loaded into a moving optical lattice 
using the truncated Wigner method for BEC dynam- 
ics 0, 0, H, El El El EH El that incorporates some 
effects of quantum fluctuations, and allows the forma- 
tion of a non-condensed fraction. Using this model we 
show that, in certain regimes of condensate dynamics, the 
anomalous heating of the condensate does not prevent its 
localization at the band edges. However, the quantum 
noise that triggers the instability will lead to different 
dynamics in every experiment. We analyse signatures of 
localisation both, in a single-trajectory realization of the 
quantum field, which is fairly similar to a "single-shot" 
experiment, and via ensemble averages. 
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II. FORMALISM 

The second-quantized Hamiltonian for a Bose gas of in- 
teracting atoms in an external trapping potential is given 
by 

H = J dx&H if + ^ J dx&&M, (1) 



Hn 



2m 



V 2 + V(x,t), 



where ^(x) is the atomic field operator that annihi- 
lates a particle at position x, m is the atomic mass and 
U = AmY?a s jm is the interaction strength, where a s is the 
s-wave scattering length. The trapping potential V(x, t) 
generally depends on both time and position. Given the 
Hamiltonian Eq. QJ, the evolution of the system's den- 
sity operator p obeys the von-Neumann equation 



dp 
~dt 



:\P,H}. 



(2) 



In general it is not possible to find analytic solutions for 
Eq. J2J), and numerical methods must be used instead. 

One manner in which to proceed is to use a phase-space 
method to represent the density operator by an expansion 
on a suitable basis, and then derive equations of motions 
for the phase space variables . In principle t he p ositive- 
P representation 01 an d its extensions [TH Il6l Il7j are 
exact methods for quantum dynamics and have had sev- 
eral successes in simulating Bose gas systems. However, 
as these methods are stochastic in nature, in general they 
suffer from an exponentially increasing sampling error 
that can limit the useful simulation time. 

In this paper we make use of another phase space tech- 
nique known as the truncated Wigner approximation. 
This was used in the context of squeezing of solitons in 
optical fibres by Drummond and Carter , and was first 
applied to Bose gases by Steel et al. [8j. We proceed 
from the Wigner function W(a(x), a*(x)) for the system 
which is derived from the density operator p(cc)_in the 
single mode case via the characteristic function 



W(a,a* 



o 0* a -0 a * Tr { e f3~J-f3*-*p} d 2p. 



(3) 



The multimode functional generalization is straightfor- 
ward, and it is possible to derive an evolution equation 
for W(a(x), a*(x)) in the form 



dW(a(x),a*(x)) 



dt 

dx — 



(4) 



h \ 5a(x) 
4 5a(x) 2 5a*(x) 



U 



Hq a(x) — U 
- c.c.l W (a(x),a*(x)) . 



If the terms involving third order functional derivatives 
with respect to a are omitted Eq. takes the form of 



a Fokker-Planck equation [6j, which can then be solved 
with stochastic methods. This approach has become 
known as the truncated Wigner approximation (TWA). 
The omission of the higher order terms in Eq. Q is justi- 
fied when there are more particles than modes in a calcu- 
lation and the simulated times are short 0, 0, IT^l Hfij . 

With the higher order derivatives neglected, the equa- 
tion of motion for the stochastic wave function a{x) that 
is equivalent to Eq. JJJ, has the same form as the usual 
Gross-Pitacvskii equation (GPE) 



ih 



da(x) 



2m 



V 2 + V + U\a(x)\ 2 ) a{x). (5) 



The difference from the Gross-Pitaevskii theory arises in 
the stochastic initial conditions, where quantum fluctua- 
tions enter the description of the initial physical state. 

In order to obtain quantum correlations of the atomic 
field Eq. ijSJ has to be solved for a large number of tra- 
jectories whose initial states a(x,t = 0) are assigned ac- 
cording to the Wigner distribution of the assumed quan- 
tum state of the BEC. In this paper we denote quan- 
tum ensemble operator averages by angle brackets e.g. 
{iff'(x)if?(x)), and averages over trajectories using over- 
lines e.g. a*(x)a(x). In the Wigner representation tra- 
jectory averages give symmetrically ordered operator av- 
erages @, for example 



a*(x,t)a(x,t) = -(&(x,t)$(x,t) + #(x, t)&(x, t)}. 

(6) 

In our numerical implementation of the Wigner method, 
we extract the condensate density and total (condensed 
and thermal) atomic density by finding the mean-field 



n (x) = | "00 (a:) I 2 , 

ritotOc) = n (x) + n t herm(x) 

= (*t(a.)*(a!)) = a *(x)a(x)- 



(7) 
(8) 
(9) 



where h x is the spacing of the computational grid |19j . 
The number of condensed or thermal atoms is therefore 
found as: 



^cond = / n (a;) da;, 



(10) 



■^Vtherm = / "therm (x) dx. 



III. 



SIMULATION PARAMETERS FOR 
SOLITON FORMATION 



We consider the quantum dynamics of a strongly 
anisotropic cigar-shaped BEC cloud loaded into a one- 
dimensional optical lattice. Provided that the conden- 
sate is tightly confined in the plane transverse to the 
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direction of the lattice, this problem can be considered 
to be effectively one-dimensional. The dimensionality of 
the equations of motion can be reduced by applying stan- 
dard procedures |2jj, which results in the following ID 
model: 

daix) fid 2 . . , ,.,\ , . , , 

l -^= [-^ + VoL(x,t)+ 1 \a(x)\ 2 'ja(x) 7 (11) 

where 7 is the rescaled one dimensional interaction 
strength 7 = 2a s /a± and length, time and energy are 
expressed in units: a± = (h/mu>±)^ , luJ 1 and huj± re- 
spectively The confining potential of the time dependent 
optical lattice is defined as Vol(x, t) = Vq sm 2 (Trx/d—vt), 
where Vq is the depth, d = A/2 is the period and v is the 
velocity of the moving lattice due to non-zero frequency 
detuning between the lattice-forming laser beams. 

Throughout this paper we give physical units of length 
and time that correspond to a condensate of 87 Rb atoms 
with atomic mass m = 1.44 x 10 _25 kg confined with the 
transverse trapping frequency uj± — 2n x 100Hz, which 
gives the unit length of a± w 1.08 /zm. We also set ir/d = 
1 which corresponds to an inter-well lattice spacing of 
~ 3.4 //m. 

We investigate the dynamics of soliton formation by 
solving Eq. (|llfl using an adaptive Runge-Kutta-Fehlberg 
method within the high level programming language 
XMDS [2lJ . We perform reliable simulations of the one- 
dimensional condensate dynamics with 8192 grid points 
on an x-range of 2ir x 512 lattice sites. This corresponds 
to a grid spacing of 0.393, approximately three times 
smaller than the condensate healing length after the soli- 
tons are formed. For our stochastic simulations we typ- 
ically use 1000 trajectories and have ascertained small 
sampling errors (less than 5%) in all ensemble averages. 



A. Initial state 

Present experimental techniques used to demonstrate 
gap solitons in a ID optical lattice Q allow access to 
the first Bloch-wave gap from the top edge of the ground 
band, where the dispersion of the BEC wavepacket is neg- 
ative but small. Therefore the strength of nonlinearity, 
proportional to 7|a| 2 , required to balance the dispersion 
and form a single soliton, is also small 3] . Hence in the 
experiment only very few (~ 350) atoms were confined 
in a single soliton. In numerical mean-field studies of a 
one dimensional system j'22f and the formation of a sin- 
gle gap soliton near the band edge with ~ 100 atoms 
was shown. This result was confirmed for BECs in the 
presence of low-energy excitations in the Bogoliubov ap- 
proximation |23j |. 

In our analysis of a condensate with a large number of 
atoms at the band edge we aim to simulate simultaneous 
generation of several gap solitons in a train, thus our ini- 
tial BEC cloud contains a few thousands of atoms if the 
natural scattering length of 87 Rb is used p|. However, 



the truncated Wigner method requires that we add a vac- 
uum noise contribution of half a particle per mode 0. 
For a typical effective computational grid size of 4096 
points this gives 2048 virtual particles. To be confident 
of the validity of the truncated Wigner method for a rea- 
sonable period of time (such that the third order deriva- 
tives in the Fokkcr-Planck equation are relatively unim- 
portant), we wish to simulate a moderately large number 
of particles so that the coherent dynamics is not over- 
whelmed by the quantum noise. 

In order to enter a regime where the TWA can be 
successfully employed in our ID simulations and soli- 
ton trains dynamically form, in addition to having a sig- 
nificant number of particles, neither the effective non- 
linearity nor the virtual population from vacuum noise 
can be too large. Only the condition on the nonlin- 
earity is physical, the others are purely practical con- 
straints arising from the very fine grids required for 
the lattice and the limitations of the TWA. We have 
found that these requirements can be reasonably met by 
reducing the interatomic interaction coupling strength 
to 7 — 10 -4 , corresponding to the scattering length 
a s = 5.4 x 10~ n m, which could in principle be achieved 
via Feshbach resonances. Wc model the initial conden- 
sate wavefunction if>o(x) as the Gross-Pitaevskii ground 
state for 55204 atoms confined in a harmonic trap of fre- 
quency uj x = 47r x 10~ 4 Hz, giving a Thomas- Fermi radius 
of about 640 /im. The initial state of our simulations is 
found by solving Eq. i|ll[l with the imaginary time prop- 
agation method in XMDS [2lj . 

We model the initial state of the Bosc g coher- 
ent state, which is realized in the Wigner distribution by 
using the initial state 

a{x,t = 0) = Mx) + -^=r ] (x), (12) 

where tpo(x) denotes the ground state of the condensate 
and r)(x) is a vacuum noise constructed from Gaussian 
random variatcs that fulfill the conditions: (i](x)r](x')) — 
and (t)*(x)t](x')) = 5(x — x'). 

For the individual simulation trajectories the initial 
state is expressed as 

a(k,t = 0)=Mk) + ^ri(k), (13) 

9 k = 6(\k\ - K max /2), 

where a denotes the spatial Fourier transform of a, i^ max 
is the maximum wave number representable on our nu- 
merical grid with 8192 points, and 6(k) is a step func- 
tion. We note that vacuum noise is only added to half 
of the momentum modes |fc| < K max /2. In our numer- 
ical simulations we ensure that significant contributions 
to the condensate evolution arises only from modes with 
\k\ < -Kmax/2, as modes with |fc| > if max /2 are poten- 
tially affected by aliasing effects due to the finite extent 
of the grid. This procedure allows us to avoid the neces- 
sity of introducing a projection method for the nonlinear 
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term into our numerical algorithm, e.g. see [3E3j an d 
references within. We have ensured in our calculations 
that initially empty high-momentum modes never have 
more that one tenth of the population of any other mode 
so that this method is justified. 



B. Loading the condensate onto the edge of the 
Brillouin zone 

In the mean-field approach, stationary states of a 
BEC are described by solutions of the classical Gross- 
Pitaevskii equation formally identical to Eq. Ijl 1|) . which 
take the form 

ipo(x, t) = ip(x) exp(-iyut), 

where /x is the chemical potential. Stationary states of 
the non-interacting BEC loaded into a periodic potential 
can be expressed as ip(x) = <fi q (x) exp(ikx), where the 
wave number k, which in the lattice rest-frame is equiv- 
alent to quasimomentum q, belongs to the first Brillouin 
zone (BZ) of the ID lattice, and <j> q (x) = cf> q (x + d) is 
a periodic (Bloch) function with the periodicity of the 
lattice. The spectrum fj,(q) of matter- waves in an op- 
tical lattice is characterised by a well-known band-gap 
diagram as shown in Fig.^a). Due to its periodicity in 
q the spectrum can be reduced to the first BZ which (in 
our units) extends over — 1 < q < 1. 

In order to access the regime of negative effective 
dispersion the condensate needs to be prepared at the 
edge of the first Brillouin zone indicated by point (2) 
in Fig. ^a). In our numerical simulations we imple- 
ment the nonadiabatic loading technique proposed and 
experimentally demonstrated by Mellish et al. [2^|. This 
method relies on the assumption that only the lowest en- 
ergy eigenstates of the lattice are involved in the loading 
process. 

In the plane wave approximation the q = and q = 1 
eigenstates in the rest frame of the lattice are: 4> q= o i x ) = 
( e - ikx + e lkx )/V2 and <f>g=i(x) = ( e - lkx - e lkx )e llJJRt /V2, 
where ton is the Rabi frequency proportional to the en- 
ergy difference between the states, </> g =o(x) is a state from 
the lower edge of the band (point (1) in Fig. Ufa)) and 
4> q =i{x) from the upper edge (point (2) in Fig.^a)). If a 
homogeneous condensate is suddenly loaded into an op- 
tical lattice moving with a positive velocity, v > 0, it 
can initially be represented in the lattice rest-frame as 
a superposition of lattice eigenstates with different en- 
ergies: tp (x) = e~ lkx = (<f> q=0 (x) + <f) q =i(x))/V2, and 
therefore will undergo Rabi oscillations between the e lkx 
and e~ lkx momentum states. After the loading, the con- 
densate can be transfered to the (f> q =i state with a sud- 
den phase shift by 8 = —tt/2 of the lattice potential, 
Vol(x) = [1 — cos(7rx/d — 2vt + 6)}, applied at the 
time t g = 3tt/(2oj r ) ~ 4.7 (7.48ms). The applied lat- 
tice displacement is equivalent to shifting the plane wave 
states by 6/2 




FIG. 1: (a) Schematics of the location (in momentum space) 
of the condensate wave packet relative to the lattice band-gap 
structure. Points (1) and (2) correspond to the modulation- 
ally stable and unstable Bloch states, (b) Rabi oscillations of 
the relative population of two k = — 1 (solid line) and k = +1 
(dashed line) momentum components of a BEC nonadiabat- 
ically loaded into an optical lattice with the height Vo = 2 
moving with the velocity v = 1, shown in the lattice rest- 
frame. The oscillations are halted by a sudden shift of the 
optical lattice by 8 — —tt/2 applied at time tg ~ 4.7 (7.48ms), 
see text. 

This reasoning can be extended to the the case of an 
inhomogeneous initial state of the BEC. As a result the 
Rabi oscillations are arrested as shown in Fig.^b). The 
initial state in our simulations is suddenly released from 
the harmonic confinement into a periodic potential. After 
the nonadiabatic loading into an optical lattice with the 
height Vq = 2, moving with the velocity v = 1, the con- 
densate undergoes Rabi oscillations, see Fig. ^b). Fol- 
lowing the 8 = —tt/2 phase shift applied to the lattice, 
the disruption of the oscillations in momentum space in- 
dicates the completion of the transfer of the condensate 
to the q = 1 momentum state, see Fig. QIb). After- 
wards, the condensate is evolved for a variable time in 
the potential of the lattice moving with a velocity v = 1 
corresponding to the Brillouin zone edge. 



IV. FORMATION OF GAP SOLITON TRAINS 

A. Localization signatures in the BEC density 

Within a relatively short time interval, the dynami- 
cal instability at the edge of the Brillouin zone starts to 
dominate the dynamics of the matter-wave and eventu- 
ally leads to condensate localization in the form of ar- 
rays of bright (gap) solitons. The development of the 
solitons from the bulk density is a complex process that 
progresses over an extended period of time. Similarly to 
the mean-field approach we find the development of 
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localised soliton-likc structures in the condensate density 
within the single trajectory TWA. As the creation of gap 
solitons requires relatively low nonlincaritics. we take an 
initial wave packet with peak density 70 (5.6 x 10 19 m~ 3 ). 

For the conditions chosen in our simulations, at the 
time corresponding to 373.3 time units (594.13ms) of evo- 
lution at the BZ edge, several structures with the char- 
acteristic beating pattern of gap solitons emerge in the 
central part of the cloud. The number of solitons devel- 
oped depends on the number of atoms in the initial BEC 
cloud and on the stochastic initial conditions. For the pa- 
rameters chosen in our single-trajectory simulations, the 
soliton train contains initially seven soliton-like peaks, 
see Fig.|2Ia). Five solitons in a sequence are separated by 
a high density region from two other solitonic structures. 
The high-density background encompassed by the train 
together with another situated to the right of the train 
do give rise to further localised structures emerging at 
later evolution times. The solitons are clearly seen in the 
density profile |a(a;)| 2 as shown in Fig.^fb), and become 
most pronounced around the time t = 511 (813.28ms). 
A single soliton from the train extends over 5 — 7 lattice 
wells. The centres of two neighbouring solitons in the 
train are separated by ~ 12 lattice periods. An average 
soliton contains about 2.5 x 10 3 atoms. 

The atoms remaining in the background accumulate 
in the vicinity of the soliton array and exhibit irregular 
modulations of density, see Fig.^a). A small fraction of 
remaining atoms is scattered in the direction opposite to 
the movement of the lattice. 
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Once formed the solitonic structures are stable for a 
long time of the simulated condensate dynamics. In the 
timespan of our simulations performed for as long as 700 
time units (1.11s), up to 14 localized soliton-like struc- 
tures form. However at longer evolution times it becomes 
more difficult to clearly distinguish localised structures 
from the background density as some of the solitons are 
characterised by very small amplitudes. The solitons de- 
velop "tails" which overlap with the background density. 

The single-trajectory simulations described above are 
useful because they provide an indication of the signa- 
tures of the localization that may be seen in a "single- 
shot" realization of an experiment [T2I Il3j . However, 
in order to determine populations of condensed and un- 
condensed atomic fractions, the stochastic formalism re- 
quires averaging over many trajectories. Due to the 
uncertainties in the position of localized peaks intro- 
duced by the noise, the averaging of the stochastic wave 
function over the ensemble does eventually lead to the 
disappeara nce of the localized structures in the density 
a*(x)a(x). Since we are interested in the enhanced ther- 
malisation of atoms that may accompany gap soliton 
train formation, we need to find another signature of the 
condensate localization that survives in the quantum en- 
semble average. 



FIG. 2: (a) Density |a(a;)| and (b) close-up of the localized 
structures in position space for a single realisation of the atom 
field, (c-f) First order density correlations in position space, 
sec Ecis. H5ll4H : (c,c) Single trajectory correlation function 
K(Ax) and (d,f) integrated second order correlation func- 
tion of the atom field averaged over 1000 trajectories G(Ax). 
Snapshots (c,d) are taken at the time t — 266 (423.35ms). 
Snapshots (a,b,e,f) are taken at the time t = 406 (646.17ms), 
after t = 401.3 (638.69ms) of evolution at the BZ edge. 



B. Density correlations 

A form of observable that can provide information 
about localization are density correlation functions. We 
define the integrated second order correlation function as 

/oo 
dx (¥(x)¥(x + Ax)$r(x)$r(x + Ax)). 
-OO 

(14) 

Since we are in the classical field regime and can with 
some validity interpret individual trajectories as corre- 
sponding to individual experiments, we can also define 
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the single trajectory correlation 

/oo 
dx \a(x + Ax)\ 2 \a{x)\ 2 , (15) 
-oo 

which fails to rigorously account for the vacuum occupa- 
tions. However, as these are random and small we expect 
coherent features to dominate. 

In Fig. |2Je) we show the shape of the correlation func- 
tion for a single trajectory, at the time after the gap 
soliton train has developed. The second order density 
correlation functions K(Ax) and G(Ax) have been nor- 
malised such that their peak densities are equal to one. 
They also exhibit very fast oscillations that are not re- 
solved on the spatial scale of the figure. 

The correlation function K(Ax) in Fig.^e) has a char- 
acteristic central peak which develops only after the soli- 
ton array has formed. The half-width of the central peak 
is about the width of a single gap soliton (roughly Ax ~ 5 
lattice periods). The development of the central peak is 
accompanied by the overall decay of the density correla- 
tions at larger Ax. Moreover, we observe clear regular 
collapses and revivals of the correlations represented by 
the sequence of maxima and minima in the envelope of 
K(Ax), as seen in Fig. Ete) . The position of the first re- 
vival {\Ax\ ~ 38) coincides with the average separation 
distance between two neighbouring gap solitons in the 
train. 

The point we stress here is that the characteristic sig- 
natures of the soliton train captured by the second or- 
der correlation function K(Ax) are also present in the 
ensemble averaging, as seen in Fig. Hff). Although the 
exact position of the solitons in every "single-shot" real- 
isation of the atom field is uncertain due to the quantum 
fluctuations, matter-wave localization occurs every time, 
because quantum noise always triggers dynamical insta- 
bility. The localization manifests itself in formation of 
soliton- like structures in the density |a(a;)| 2 . The main 
advantage of the density correlation function is that it 
captures the average size of these localized structures 
(and spacing between them) instead of the exact loca- 
tion of the solitons. Hence the density correlations are 
less sensitive to soliton position fluctuations between dif- 
ferent trajectories than the averaged density. 

Similarly to the single trajectory correlation function, 
G(Ax) also acquires the characteristic peak values when 
the localization occurs and the probability of finding two 
atoms separated by a certain distance in the lattice is the 
highest. This happens either when Ax is equal to zero 
or to the separation distance between centres of neigh- 
bouring solitons (Ax as 38). Between the central and a 
neighbouring maximum of G(Ax), a local minimum of 
the correlation function at |Ax| ~ 18, which correspond 
to "dips" of the matter- wave density between two neigh- 
bouring solitons, can be seen in Fig.^f). 
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FIG. 3: (a) Density profile a*(x)a(x) and (b,c) second or- 
der correlation functions G(Ax), averaged over 1000 trajecto- 
ries, for the case of initial phase imprinting. Snapshots (a,c) 
are taken at t — 119 (189.39ms), which is after t = 114.3 
(181.91ms) of evolution at the BZ edge. Snapshot of the cor- 
relation function before the gap soliton train develops (b) is 
taken at t = 56 (89.13ms). 



V. PHASE IMPRINTING 

If the dynamical instability is seeded only by the quan- 
tum noise of the initial state, like in Eq. I|12[l. the posi- 
tions of the localised structures at each realization of the 
experiment would be random. If one would like to gener- 
ate and manipulate the emerging solitons in a controlled 
fashion, it would be helpful to accurately forecast their 
position. This can be done by the imprinting of a periodic 
phase onto the initial BEC cloud before the condensate 
is loaded into a moving lattice. In this section, we will 
show that the merits of the periodic imprinting of an ini- 
tial phase onto the condensate are two- fold. First, we 
show that the phase-imprinting leads to periodic density 
modulations of the condensate which facilitates the de- 
velopment of modulational (dynamical) instability and 
significantly reduces the time scale of the soliton train 
formation. Secondly, seeding of the periodic phase (den- 
sity) modulations leads to the development of solitons at 
fixed locations. 

The phase imprinting techniques, well developed for 
BECs [221 128| , were also used to place low-atom number 
condensates at the edge of the first Brillouin zone [23l| . 
They result in an engineered spatially dependent phase 
factor exp(icj>(x)) of the BEC cloud. Here wc analyse the 
signatures of the localization observables if the phase of 
the initial BEC is regulary modulated as <f>(x) = cos(S-x), 
where we take 6 = 0.1 for the results we plot. Our sim- 
ulation sequence is as follows: First, the initial state in 
the harmonic trap is obtained for the same parameters as 
previously. Next, the periodic modulation of the phase 
is imprinted onto the BEC cloud. Then, similarly to the 
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previously described scheme, the real time evolution be- 
gins, and the condensate is nonadiabatically loaded into 
the moving optical lattice at the Brillouin zone edge. Fi- 
nally, the condensate is left to evolve in the moving lattice 
for a variable time up to t = 700 (1.11s), which is equiv- 
alent to t = 695.3 (1.1s) of evolution at the BZ edge. 

As a result of the phase imprinting, strong den- 
sity modulations as well as striking regularities within 
the correlation function develop at very early stages of 
the evolution. Periodic density modulation usually at- 
tributed to triggering modulational (dynamical) insta- 
bilities can be clearly seen already at the time t = 14 
(22.28ms). Because of these, the process of the gap soli- 
ton array formation is greatly accelerated. The solitons 
can be clearly visible in a single-trajectory simulation af- 
ter only t = 93.3 (148.49ms) of evolution at the BZ edge. 
There is also a larger number of solitons in a train (ini- 
tially there are as many as 18 localised structures). The 
number of emerging solitons and their separation does 
not depend on the period of phase modulation S as long 
as the wavelength A = 2tt/S corresponding to the phase 
perturbation is larger than the characteristic scale for gap 
solitons (soliton width). 

The soliton-like structures are stable up to t = 175 
(278.52ms). Around t ~ 200 their dynamics become 
more complicated and due to their interactions the num- 
ber of solitons in the array changes [2|j. At even longer 
evolution times it becomes more difficult to clearly iden- 
tify localised structures from the background, however 
about 14 soliton-like structures survive until the end of 
"single-shot" simulations corresponding to 1.1s. 

Most remarkably, the formation of the solitons at fixed 
positions manifests itself in the visibility of the soliton 
train in the condensate density profile even after ensem- 
ble averaging as shown in Fig. |3Ja). The robustness of 
the position of the matter- wave localisation achieved with 
the aid of phase imprinting reveals itself not only in the 
fact that the solitons survive in the quantum ensemble 
average, but also in the much more pronounced peri- 
odic structure of the second order correlation function, 
as shown in Fig. |3fc) for the case of 1000 trajectories. 
Clear periodic revival of the second order (density) cor- 
relations can be seen. 



VI. COMPARISON WITH THE 
GROSS-PITAEVSKII MODEL 

It has been shown [j| that the nonlinear localisation of 
matter-waves due to dynamical instability at the edge of 
the first Brillouin zone and the development of gap soli- 
ton trains is a nonlinear effect which can be described by 
the Gross-Pitaevskii model, i.e. Eq. (JSJ without the ad- 
dition of noise to the initial state. Those studies were 
performed for a condensate adiabatically loaded onto 
the edge of the first Brillouin zone, therefore the mod- 
eled dynamics was extended over a long period of time. 
Nonetheless even a comparison of our stochastic simu- 



lations with the mean-field evolution within the much 
faster scheme of nonadiabatic loading shows, that for the 
same parameters the localization arises about four times 
earlier in the quantum theory. We note that if initial 
phase imprinting is used, it governs the time-scale for 
the instability and in this scenario the train appears at 
similar times with and without the quantum noise. We 
stress, that the results presented here are more accurate 
than simulations of the Gross-Pitaevskii equation, which 
rely on the exponential growth of numerical inaccuracies 
to show the physical effects that we have demonstrated 
to be in fact seeded by irreducible quantum noise. 



VII. VALIDITY OF SIMULATIONS 

Unfortunately it is difficult to estimate the simula- 
tion time for which the truncated Wigner approxima- 
tion is valid and the numerical results are quantita- 
tively reliable. However, for the parameters we have 
chosen the requirement of modeling significantly more 
particles (55000) than effective modes (4096) is fulfilled. 
These numbers put our calculations into the classical 
field regime, where the quantum noise is effectively a 
seed for more complicated, chaotic dynamics arising from 
the dynamical instability ^lj. These result in heating 
and a non-zero temperature final state, which can not be 
treated by the mean-field model, see section IVIIII Once 
begun, the dynamics are dominated by nonlinear inter- 
actions between highly occupied modes so that even if 
our method gives results that are quantitatively inaccu- 
rate beyond a certain time, there will be a much longer 
time-scale for which they are qualitatively correct. We 
certainly expect that the formation of soliton trains that 
we have demonstrated will persist in any experiment that 
can be peformed given our initial starting point. 

VIII. ANOMALOUS HEATING 

As described earlier, in addition to gap soliton train 
formation a condensate evolving at the ed ge o f the first 
BZ exhibits enhanced loss of atoms pH l30j| . Previ- 
ously 0,H3,|^|, the onset of the dynamical instability of 
the Bloch state at the BZ edge was conclusively linked 
to the growth of the thermal fraction. We stress, that 
anomalous heating can also be observed for condensates 
with quasimomenta far from the BZ edge |30j |. however, 
this process is triggered by the energetic (Landau) in- 
stabilities that are initiated by the presence of a finite 
thermal component |30j and facilitate decay of the con- 
densate into the lowest energy state. In contrast, in our 
simulations the initial state is prepared as a pure coher- 
ent state BEC. 

For the parameter regime used in our simulations, the 
formation of gap soliton trains is accompanied by the en- 
hanced thermalisation of atoms in the cloud. In Fig. 0] 
we present the transfer of atoms from the condensed to 



8 



x 10 




FIG. 4: Loss of BEC atoms (a) and growth of the thermal 
fraction (b) during t — 4.7 (7.48ms) of nonadiabatic loading 
followed by t = 695.3 (1.1s) of the evolution at the BZ edge. 
Shown are numbers of condensed and uncondensed atoms dur- 
ing the evolution without (solid lines) and with (dashed lines) 
the initial phase imprinting. 
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FIG. 5: (a,b) Close-ups of averaged density profiles 
a*(x)a(x) for the scheme with addition of the phase imprint- 
ing. Shown are snapshots from two different time samples 
indicated by matching red circles in (c). BEC profiles corre- 
spond to: (a) t = 98 (155.97ms) or 93.3 (148.49ms) of the 
evolution at the BZ edge, and (b) t = 126 (200.54ms) or 
t = 121.3 (193.05ms) of the evolution at the BZ edge, (c) 
Same as in Fig.0]but for shorter evolution time, up to t — 300 
(477.46ms). 



the uncondensed fraction of the Bose gas during the time 
evolution within our scheme. Also shown is the thermal- 
isation for the case of phase imprinting applied initially 
onto a(x) for every trajectory, see Fig. 01 and Fig.0 

The phase modulation results in a smaller loss of atoms 
from the condensate at the long time scales of our sim- 
ulations, see Fig. 0] However initially at the times when 



solitons form t ~ 120, it leads to an increased rate of 
atom loss, see Fig. [5J which is consistent with the obser- 
vation that the phase modulation accelerates the dynam- 
ical instability. It can be seen in Fig. [31 that when the 
excitations of the condensate increase, in the case with 
initial phase modulations, the bulk density collapses and 
localisation commences. In Fig. EJb), showing the time 
sample corresponding to t = 126 (200.54ms), the devel- 
opment of the soliton arrays can be clearly seen even in 
the averaged density profile. 



IX. PARAMETER REGIMES FOR THE 
LOCALISATION 

Furthermore, we have investigated the localisation in 
different parameter regimes. To this end we have mon- 
itored the condensate dynamics for higher interaction 
strengths: with coupling 7 increased to S7, keeping all 
other parameters fixed. This reduces the peak-density at 
t = by a factor of s (\ij} (t = 0,x)\ 2 = /zs"^ 1 ). We 
found that for s = 10 wc do still observe the localisa- 
tion of the matter-wave but due to the amount of noise 
present the pattern is much less regular than in the case 
presented in Fig. [2J We suspect this parameter regime 
to be the borderline of the applicability of the Wigner 
method. We draw such a conclusion from probing the 
BEC dynamics for s = 100. In this situation we observe 
that any signatures of the localisation are overwhelmed 
by the noise. 

In addition, we have investigated the dynamics with 
m times higher nonlincarities but fixed densities: The 
coupling strength 7 and the chemical potential fi were 
increased by a factor of to and adequately the trapping 
frequency of the HO was increased by y/m, setting the 
Thomas-Fermi radius constant. In the case when to = 10 
our single-trajectory simulations do show a fragmenta- 
tion of the matter- wave: We observe the development of 
high density spikes which we attribute to the anomalous 
heating effect in the "single-shot" realisation. Indeed in 
the quantum ensemble average the heating becomes more 
dominant. It shows a much greater amount of thcrmalisa- 
tion at early stages of the evolution (t <C 100). However, 
in the case m = 100 thermalisation is so intensive, that 
no mark of fragmentation of the matter-wave is visible in 
a single-trajectory simulation. 

For to = 10 the localisation signatures in the corre- 
lation functions (K(Ax) and G(Ax)) become very weak 
but the minima in the envelope arc preserved. On the 
other hand for to = 100 all the signatures completely 
disappear. 



X. CONCLUSIONS 

Within the phase-space truncated Wigner method we 
have investigated the dynamics of matter-waves in a mov- 
ing lattice beyond the onset of dynamical instability. We 
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have demonstrated localisation of an atomic field within 
the single and multitrajectory treatment of the TWA. 
We have shown that a train of strongly localised matter- 
wave gap solitons can be generated as a result of dynam- 
ical instability of a BEC nonadiabatically loaded into a 
moving optical lattice. These instabilities are triggered 
by quantum noise. The localization happens despite the 
increased thcrmalisation of the condensate at the edge of 
the Brillouin zone. 

We proposed the acceleration of gap soliton array for- 
mation from moderate-size atomic clouds at the edge of 
the band by a phase imprinting technique. In addition we 
have demonstrated that the initial phase imprinting onto 
the BEC cloud leads to the localisation of matter-waves 
at fixed positions. This method ensures the formation of 



spatially regular soliton trains at short evolution times. 
Density-density correlation functions are shown to dis- 
play clear signatures of the condensate localization. 
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